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ABSTRACT 

In the context of difference image analysis (DIA), we present a new method for determin- 
ing the convolution kernel matching a pair of images of the same field. Unlike the standard 
DIA technique which involves modelling the kernel as a linear combination of basis functions, 
we consider the kernel as a discrete pixel array and solve for the kernel pixel values directly 
using linear least-squares. The removal of basis functions from the kernel model is advanta- 
geous for a number of compelling reasons. Firstly, it removes the need for the user to specify 
such functions, which makes for a much simpler user application and avoids the risk of an in- 
appropriate choice. Secondly, basis functions are constructed around the origin of the kernel 
coordinate system, which requires that the two images are perfectly aligned for an optimal 
result. The pixel kernel model is sufficiently flexible to correct for image misalignments, and 
in the case of a simple translation between images, image resampling becomes unnecessary. 
Our new algorithm can be extended to spatially varying kernels by solving for individual pixel 
kernels in a grid of image sub-regions and interpolating the solutions to obtain the kernel at 
any one pixel. 
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1 INTRODUCTION 

Difference image analysis (DIA) has rapidly moved to the forefront 
. , of modem techniques for making time-series photometric measure- 
. ments on digital images. The method attempts to match one image 
' to another by deriving a convolution kernel describing the changes 
5^ i m m e point spread function (PSF) between images. When applied 
to a time-series of images using a high signal-to-noise reference 
image, the differential photometry that can be performed on the 
difference images regularly provides superior accuracy to more tra- 
ditional profile-fitting photometry, achieving errors close to the the- 
oretical Poisson limits. Moreover, DIA is the only reliable way to 
analyse the most crowded stellar fields. 

One will find DIA in use in many projects st udying object vari- 
ability. For example , microlensing searches (e.g. Bon d et alj|200ll ; 
Wozniak et alj200lh have been revolutionised by the ability of DIA 
to deal with exceptionally crowded fields, and surveys for t ransiting 
planets (e.g. iBramich et alj|2005l ; iMocheiska et alj|2005h looking 
for small ~1% photometric eclipses have benefited substantially 
from the extra accuracy obtained with this method. Also, DIA is 
not limited to stellar photometry as illustrated by the discovery of 
light e choes from three ancient supernovae in the Large Magellanic 
Cloud teest et alj|2005t) . 

The first at tempts at image subtrac tion are summarised in the 
introduction of lAlard & Lup ton ( 1998) (from now on AL98) and 



are based on trying to determine the convolution kernel by taking 
the ratio of the F ourier transforms of matc hing bright isolated stars 
on each image dTomanev & Crottj Il996h . Development of DIA 
reached an important landmark in AL98 with their algorithm to 
determine the convolution kernel directly in image space (rather 
than Fourier space) from all pixels in the images by decompos- 
ing the kernel onto a set of basis functions. The algorithm is very 
successful and efficient, and with the exten sion to a space- varying 
kernel solution described in lAlard feOOCh (from now on AL00), 
the method has become the current standard in DIA. In fact, all 
DIA packages use the associated softw are package ISIS2. J] (e.g. 
Wozniak 2000; Gossl & Riffeser 2002), or are implementations of 
the Alard algorithm (e.g. lBond et alJuOOll) . We refer to the method 
described in AL98 and AL00 as the Alard algorithm. 

In this letter we suggest a change to the main algorithm to de- 
termine the convolution kernel that retains the linearity of the least- 
squares problem and yet is simpler to implement, has fewer input 
parameters and is in general more robust (Section 2). We compare 
our algorithm directly to the Alard algorithm (Section 3), and sug- 
gest more techniques that increase the quality of the subtracted im- 
ages. We conclude in Section 4. 
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2 A NEW APPROACH TO THE KERNEL SOLUTION 
2.1 Motivation 

Consider a pair of registered images of the same dimensions, one 
being the reference image with pixels Rij, and the other the cur- 
rent image to be analysed with pixels /y, where i and j are pixel 
indices refering to the column i and row j of the image. Ideally the 
reference image will be the better seeing image of the two and have 
a very high signal-to-noise ratio. This can be achieved in practice 
by stacking a set of best-seeing images. 

As with the method of AL98, we use the model 



Mij = [R®K]i 



Bi 



(1) 



to represent the current image Ty, where we wish to find a suit- 
able convolution kernel K and differential background Bij. For- 
mulating this as a least-squares problem, we want to minimise the 
chi-squared 



X 



Mi. 



(2) 



where the cry represent the pixel uncertainties. 

At this point in the Alard algorithm, the problem is converted 
to standard linear least-squares by decomposing the kernel K onto 
a set of gaussian basis functions, each multiplied by polynomials 
of the kernel coordinates u and v, and by assuming that the differ- 
ential background Bij is represented by a polynomial function of 
the image coordinates x and y. Spatial variation of the convolution 
kernel is modelled by further multiplying the kernel basis functions 
by polynomials in x and y. 

This method has a major drawback in that it assumes that the 
chosen kernel decomposition is sufficiently complex so as to model 
in detail the convolution kernel. How do we know that we are mak- 
ing the correct choice of basis functions? Different situations may 
require different combinations of basis functions of varying com- 
plexity. In fact, a feature of all current DIA packages (which are 
all based on the AL98 prescription for kernel basis functions) is 
the requirement that the user defines the number of gaussian basis 
functions used, their associated sigma values and the degrees of the 
modifying polynomials. This sort of parameterisation can end up 
being confusing for the user, and require a large amount of experi- 
mentation to obtain the optimal result for a specific data set. 

2.2 Solving For A Spatially Invariant Kernel Solution 

With this motivation, we have developed a new DIA algorithm in 
which we make no assumptions about the functional form of the 
basis functions representing the kernel. Considering a spatially in- 
variant kernel, we represent the kernel as a pixel array Ki m with 
Nk pixels where I and m are pixel indices corresponding to the 
column I and row m of the kernel. We also define the differential 
background as some unknown constant Bo- Hence we may rewrite 
equation |QJ as: 



Mi, 



of x 2 with respect to each of the parameters Ki m and Bo is equal to 
zero. Performing the Nk + 1 differentiations and rewriting the set 
of linear equations in matrix form, we obtain the matrix equation 
Ua = b with: 



U v 



E* 
E* 



E« 



fl ( , +0(j+m) fi ( , +i0(jW) for x < p < Nr . and 1 < g < Nr 



°y 

^ 



for p = Nk + 1 and 1 < q < N K 
for 1 < p < N K and q = N K + 1 
for p = q = Nk + 1 



K lm for 1 < p < N K 
B forp = N K + l 



bp = 



J ij-"(i + Q(j + m) 



for 1 < p < N K 
forp = Nk + 1 



(4) 



where p and q are generalised indices for the vector of unknown 
quantities a, with associated kernel indices (l,m) and (I' ,m') re- 
spectively. Finding the solutionflfor Ki m and Bo requires the con- 
struction of the matrix U and vector b, inverting U and calculating 
a = U" 1 b. 

Every pixel on both the reference image and current image has 
the potential to be included in the calculation of U and b. However, 
we ignore bad/saturated pixels on both images, and also any pixels 
on the current image for which the calculation of the correspond- 
ing model pixel value includes a bad/saturated pixel on the refer- 
ence image. This implies that a single bad/saturated pixel on the 
reference image can discount a set of pixels equal to the kernel area 
on the current image. Hence bad/saturated pixels on the reference 
image should be kept to a minimum, and excessively large kernels 
should be avoided. 

The kernel sum P — ^ m Ki m is a measure of the mean 
scale factor between the reference image and the current image, and 
consequently it includes the effects of relative exposure time and at- 
mospheric extinction. We refer to P as the photometric scale factor. 
Although it is not essential, we suggest that a constant background 
estimate is subtracted from the reference image before solving for 
the kernel and differential background since this will minimse any 
correlation between P and Bo. 

Finally, we mention that a difference image Dy is defined as 



Di 



Mij. Assuming that most objects in the reference 



image are constant sources, then a difference image will consist of 
random noise (mainly Poisson noise from photon counting) except 
where a source has varied in brightness or the background pattern 
has varied. Sources that are brighter or dimmer at the epoch of the 
current image relative to the epoch of the reference image will show 
up as positive or negative flux residuals, respectively, on the differ- 
ence image. These areas may be measured to yield a difference flux 
for each object of interest. 



lmR(i + l)(j + m) 



(3) 



This equation has Nk + 1 unkowns for which we require a solution. 
Note that the kernel may be of any shape that includes the pixel 
Koo, and so to preserve symmetry in all directions, we adopt a 
circular kernel (instead of the standard square shape). 

In order to solve for K\ m and Bo in the least-squares sense, we 
note that the x' 2 m equation (O is at a minimum when the gradient 



Iteration is required for a self-consistent solution for K\ m and Bo since 
the solution depends on the pixel variances <r^ which in turn depend on the 



image model values Afjw. See Section l231 
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2.3 Uncertainties Arise From The Model, Not The Data 

We take the following standard CCD noise model for the pixel vari- 
ances: 

2 ol Mij 

where o"o is the CCD readout noise (ADU), G is the CCD gain 
(e~/ADU) and Fij is the master flat field image. Note that the o-ij 
depend on the image model My and consequently, fitting My be- 
comes an iterative process. Note also that we assume that the refer- 
ence image Rij and master flat field image Fy are noiseless since 
these are high S/N ratio images. Finally, if the current image was 
registered with the reference image via a geometric transformation, 
then the flat field Fij that is actually used in the noise model must 
be the result of the same transformation applied to the original mas- 
ter flat field. 

In order to calculate an initial kernel and differential back- 
ground solution, we set the My to the image values Uj. In subse- 
quent iterations, we use the current image model to set the cry as 
per equation [5] We also employ a 3a clip algorithm during the it- 
erative model fitting process in order to prevent outlier image pixel 
values from entering the solution. After each iteration, we calculate 
the absolute normalised residuals ry = |-Dy/<7y| for all pixels. 
Any pixels with ry > 3 are ignored in subsequent iterations. The 
iterations are stopped when no more image pixels are rejected and 
at least two iterations have been performed. 

2.4 Solving For A Spatially Variant Kernel Solution 

In extending our new method to solving for a spatially variant ker- 
nel solution, we preserve flexibility by splitting the image area into 
an N x > 1 by N y > 1 grid of sub-regions and solving for the ker- 
nel and differential background in each sub-region. The coarse grid 
of kernel and differential background solutions may be interpolated 
to yield the solution corresponding to any given image pixel. In this 
way we make no assumptions about how the kernel and differential 
background vary across the image area. This is in contrast to ALOO, 
whose method employs an extension of the kernel basis functions 
by further multiplication by polynomials in x and y, and therefore 
requires two more input parameters from the user, namely the de- 
grees of the polynomials describing the spatial variation of the ker- 
nel and the differential background. 



3 COMPARISONS WITH THE ALARD ALGORITHM 
3.1 Initial Tests 

To illustrate the potential advantages of our new kernel solution 
method over that of AL98, we carry out a set of simple tests on a 
1024x 1024 pixel CCD image of the globular cluster NGC1904. In 
each test we use the original image as the reference image Rij and a 
transformed version of the original image as the current image Iy, 
where the transformations employed are simple, spatially invariant 
and typical of astronomical imaging. We attempt to solve for the 
kernel using our new method, which is implemented in a software 
package called DANDIA (Bramich in prep.), and we compare the 
solution to that obtained using the ISIS2.2 software from ALOO. 
We use the ISIS2.2 default parameters specifying 3 gaussian basis 
functions of a = 0.7, 2.0, 4.0 pix with modifying polynomials of 
degree 6, 4 and 3, respectively. For both software packages, we 



choose to solve for a spatially invariant kernel of size 27 x 27 pixels, 
and a constant differential background. 

The better the match between the convolved reference im- 
age and the current image, the smaller the value of the quantity 
S 2 — Yli j D 2 j- We guage the relative quality of the kernel solu- 
tions by calculating the noise ratio Sjsis /^DANDIA where Sjsis 
and ^DANDIA are vames of S calculated for a small 80x80 pixel 
sub-region using ISIS2.2 and DANDIA, respectively. 

The results of the tests described below are shown in Figure[T] 

(i) In test A, the current image has been created by shifting the 
reference image by one pixel in each of the positive x and y spatial 
directions, without resampling. The corresponding kernel should be 
the identity kernel (central pixel value of 1 and elsewhere) shifted 
by one pixel in each of the negative u and v kernel coordinates. 
DANDIA recovers this kernel to within numerical rounding errors 
whereas ISIS2.2 recovers a peak pixel value of 0.995 with other 
absolute pixel values of up to 0.004. Consequently the residuals in 
the ISIS2.2 difference image are considerably worse than those for 
DANDIA, and the noise ratio is Sjsis/Sdandia ~ 26190. 

(ii) In test B, the current image has been created by convolv- 
ing the reference image with a gaussian of FWHM 4.0 pix. Both 
DANDIA and ISIS2.2 recover the kernel successfully, but DAN- 
DIA out-performs ISIS2. 2 with Sjsjs/SrjANDIA ~ 1510. 

(iii) In test C, we shifted the reference image by half a pixel in 
each of the positive x and y spatial directions to create the current 
image, an operation that required the resampling of the reference 
image. We used the cubic O-MOMS resampling method (see Sec- 
tion |3,2| >, ISIS2.2 fails to reproduce the highly complicated kernel 
matching the two images, whereas DANDIA does a nearly perfect 
job. The noise ratio is Sisis/SbANDIA ~ 18450. 

(iv) In test D, we simulate a telescope jump by setting 
Iy = (0.6 x Rij) + (0.4 x R'ij) where R\j is a resam- 
pled version of the reference image shifted by 3.5 pixels in each 
of the positive x and y spatial directions. The corresponding kernel 
is a combination of the identity kernel and a shifted version of the 
kernel from test C. DANDIA accurately reproduces this kernel with 
a central pixel value of 0.60015 whereas ISIS2.2 produces a poor 
approximation of the kernel with a central pixel value of 0.63 1 . The 
noise ratio is <Sisis /^DANDIA ~ 31000. 

It is evident that the gaussian basis functions used in ISIS2.2 
limit the flexibility of the kernel solution to modelling kernels that 
are centred near the kernel centre and that have scale sizes similar 
to the sigmas of the gaussians employed. It is only in test B that 
ISIS2.2 can closely model the kernel, simply because the kernel it- 
self is a gaussian. Tests A, C & D show how ISIS2.2 is unable to 
model sharp, complicated and off-centred kernels. DANDIA does 
not suffer from any of these limitations since it makes no assump- 
tion about the kernel shape, and hence it performs superbly in all 
of the above tests. 

3.2 Image Resampling 

In Section 2, we make the assumption that the reference image and 
current image are registered, which implies that one of the images 
has been transformed to the pixel coordinate system of the other 
image, usually via image resampling. Ideally one should transform 
the reference image to the current image since the reference im- 
age forms part of the model. In this way, the pixel variances in the 
current image are left uncorrelated from pixel to pixel. However, 
most implementations of DIA transform the current image to the 
coordinate system of the reference image using image resampling. 
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Figure 1. Shown are four difference imaging tests A-D. In each panel corresponding to a test, there are three pairs of images, where each pair is shown using 
the same intensity scale indicated by the graduated colour bar. In the left most pair, the reference image and current image sub-regions are shown on the left 
and right, respectively. In the middle pair, the DANDIA and ISIS2.2 difference images are shown at the top and bottom, repsectively. In the right most pair, 
the corresponding DANDIA and ISIS2.2 kernel solutions are shown at the top and bottom, respectively. 



We suggest two improvements to this methodology. Firstly, if 
resampling is to be employed, one should use an optimal resam- 
pling method. We employ the cubic O-MOMS (Optimal Maximal- 
Order-Minimal-Support) basis function for resampling, which is 
constructed from a linear combination of the cubic B-spline func- 
tion and its derivatives. The O-MOMS class of functions have the 
highest approximation or der and smalles t approximation error con- 
stant for a given support dBlu et alfeOQll) . 

Secondly, our kernel model does not use basis functions that 
are functions of the kernel pixel coordinates. Consequently, for two 
images that require only a translation to be registered, the image re- 
sampling is incorporated in the kernel solution, avoiding the prob- 
lem of correlated pixel noise. DIA is used extensively for extract- 
ing lightcurves of objects in time-series images, which usually only 
have a small pixel shift between images. By translating the current 
image to the reference image by an integer pixel shift, avoiding im- 
age resampling, the kernel solution process can do the rest of the 
job of matching the reference image to the current image. 

3.3 Final Tests 

We now test our new algorithm on a pair of 1024x 1024 pixel im- 
ages of NGC1904 from the same camera with FWHMs of ~3.2 pix 
and ~4.9 pix. Using matching star pairs, we derive a linear transfor- 
mation between the images that consists of a translation with neg- 
ligible rotation, shear and scaling. From the calibration images, we 
measure a gain of 1.48 e~/ADU and a readout noise of 4.64 ADU, 
and we construct a master flat field for use in the noise model. On 
the left of Figure[2] we present lOOx 100 pixel cutouts of the refer- 
ence image (the better seeing image) and the current image. 

When calculating the \ 2 of the difference images, we use a 



modified version of equation[5]to account for the noise contribution 
from the single-exposure reference image: 

where Kimij is the space variant kernel and / is a factor correct- 
ing for the noise distortion from resampling the reference image. 
The value of / depends on the resampling method used and the 
coordinate transformation applied. We calculate / by generating 
a 1024x 1024 pixel image of values drawn from a normal distri- 
bution with zero mean and unit sigma, resampling the image using 
the same method and transformation as that applied to the reference 
image, and then fitting a gaussian to the histogram of transformed 
pixel values, the sigma of which indicates the value of /. For cubic 
O-MOMS resampling and the transformation between our two test 
images, we obtain / = 0.884. 

Our first pair of tests involves registering the images by re- 
sampling the reference image via cubic O-MOMS and then using 
DANDIA (test E) and ISIS2.2 (test F) to generate difference im- 
ages. For DANDIA, we solve for an array of circular kernels corre- 
sponding to a 10 x 10 grid of image sub-regions, where each kernel 
contains 317 pixels. The kernel used to convolve each pixel on the 
reference image is calculated via bilinear interpolation of the array 
of kernels. The results of test E are displayed in the upper middle 
panel of Figure[2] where we show the difference image normalised 
by the pixel noise from equation [6] with a linear scale from -2 to 2. 
Two variable stars are visible (RR Lyraes) and the cosmic ray from 
the reference image has created a negative flux on the difference 
image. In the same panel we plot the histogram of normalised pixel 
values overlaid with a gaussian fit, and calculate a \ 2 ~ 9439, 
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Figure 2. Shown are four more difference imaging tests E-F. On the left we show lOOx 100 pixel cutouts from the reference image and current image. The 
remaining panels show corresponding difference images normalised by the pixel noise model (equation[6j with a linear scale from -2 to 2, and histograms of 
the normalised pixel values overlaid with a gaussian fit. 



ignoring the small pixel areas including the variable stars and the 
cosmic ray (250 pix). The lOOx 100 pixel cutout corresponds to one 
image sub-region used to determine a kernel solution and hence we 
may calculate a reduced chi-squared x 2 /N^ov = 1.00 by assum- 
ing JV D0F = 9750 - 318. 

For ISIS2.2 we solve for a spatially variant kernel of degree 2 
with a spatially variant differential background of degree 3 in addi- 
tion to the other default kernel basis functions (see Section [3~Tl 328 
free parameters). The results of test F are shown in the upper right 
panel of Figure [2] We obtain x 2 ~ 9838, and assuming ~3 free 
parameters per image sub-region, we obtain x 2 /A'dof = 1.01. 

Tests G & H involve registering the images to within 1 pixel 
by translating the reference image via an integer pixel shift. Then 
we apply DANDIA (test G) and ISIS2.2 (test H) to obtain kernel 
solutions, avoiding the use of resampling. For DANDIA we obtain 
X 2 ~ 9373, and for ISIS2.2 we obtain x 2 ~ 9739, with corre- 
sponding xV-^dof of 0.99 and 1.00, respectively (see Figure[2]l. 

Visually, the normalised difference image cutouts in Figure [2] 
are very similar, and differences are only noticeable after detailed 
scrutiny. However, the x 2 analysis reveals that our algorithm per- 
forms considerably better than the Alard algorithm (test E performs 
0.60(7 better than test F, and test G performs 0.38<r better than test 
H), and that image resampling degrades the difference images (test 
G performs 0.48er better than test E, and test H performs 0.70ct bet- 
ter than test F). The highest quality difference image was produced 
by using DANDIA on the two images aligned to within 1 pixel but 
without resampling (test G, which performs 1.08a better than test 
F). 



4 CONCLUSIONS 

We have presented a new method for determining the convolution 
kernel matching a best-seeing reference image to another image of 
the same field. The method involves modelling the kernel as a pixel 
array, avoiding the use of possibly inappropriate basis functions, 
and eliminating the need for the user to specify which basis func- 
tions to use via numerous parameters. For images that require a 



translation to be registered, the kernel pixel array incorporates the 
resampling process in the kernel solution, avoiding the need to re- 
sample images, which degrades their quality and creates correlated 
pixel noise. Kernels modelled by basis functions may only partly 
compensate for sub-pixel translations since the basis functions are 
centred at the origin of the kernel coordinates. 

We have shown that our new method can produce higher qual- 
ity difference images than ISIS2.2. Ideally the reference image 
should be aligned with the current image, preferably without re- 
sampling, but using O-MOMS resampling when necessary. The 
flexibility of our kernel model allows the construction of differ- 
ence images for telescope jumps, or trailed images, which is where 
ISIS2.2 fails. These improvements have important implications for 
time-series photometric surveys. Better quality difference images 
implies more accurate lightcurves, and the increased kernel flex- 
ibility will lead to less data loss due to telescope tracking and/or 
focus errors. 
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